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Quantum Monte Carlo and semiclassical methods are used to solve two and four site cluster 
dynamical mean field approximations to the square lattice Hubbard model at half filling and strong 
coupling. The energy, spin correlation function, phase boundary and electron spectral function 
are computed and compared to available exact results. The comparision permits a quantitative 
assessment of the ability of the different methods to capture the effects of intersite spin correlations. 
Two real space methods and one momentum space representation are investigated. One of the two 
real space methods is found to be significantly worse: in it, convergence to the correct results is 
found to be slow and, for the spectral function, nonuniform in frequency, with unphysical midgap 
states appearing. Analytical arguments are presented showing that the discrepancy arises because 
the method does not respect the pole structure of the self energy of the insulator. Of the other 
two methods, the momentum space representation is found to provide the better approximation 
to the intersite terms in the energy but neither approximation is particularly acccurate and the 
convergence of the momentum space method is not uniform. A few remarks on numerical methods 
are made. 

PACS numbers: 71.10.Fd, 75.10.-b 



I. INTRODUCTION 

"Strongly correlated" materials^ pose one of the out- 
standing challenges in condensed matter physics. These 
materials exhibit a wide range of interesting and po- 
tentially useful properties including high temperature 
superconductivity 2 - and magnetism with very high spin 
polarization;- however, in these classes of material the 
electron-electron and electron-lattice interactions are so 
strong that the conventional approach (using density 
functional theory to compute bands and then using per- 
turbative methods to treat residual interactions between 
quasiparticles) fails. Developing a reliable, material- 
specific theoretical framework for determining the behav- 
ior of strongly correlated compounds is an outstanding 
challenge to materials theory. 

The development of single-site dynamical mean field 
theory^ was a fundamental step forward in correlated 
electron science. In this approach one approximates the 
momentum and frequency dependent self energy S(p,oj) 
by a momentum-independent function S(p, — > 
This approximation allows the construction of a nonper- 
turbative and computationally tractable theoretical pro- 
cedure for computing physical properties: because it is 
a function of only one frequency variable the self energy 
may be viewed as the self energy of a single-site "quan- 
tum impurity model" , with the parameters of the model 
specified by a self consistency condition. The approach 
works very well for situations (including the Mott transi- 
tion in electronically three dimensional materials^ the 
"double exchange" physics important for colossal ma- 
gentoresistance manganites,— and the basic physics of 
heavy fermion compounds^), in which Galilean invari- 
ance is strongly broken and the dominant physics is on- 



site. However, in wide classes of interesting materials, in- 
tersite correlations play an important role in the physics. 
Examples include the high temperature superconductors, 
where the predictions of the single-site dynamical mean 
field theory have been shown to disagree strongly with 
data on the evolution with doping of quasiparticle ve- 
locity and 'Drude' optical weight^ and the orbital or- 
der/polaron glass physics of the manganites. 9 Extension 
of the dynamical mean field method to include intersite 
correlations is therefore an important issue. 

The single-site dynamical mean field theory involves 
the mapping of a lattice model onto a single-site quan- 
tum impurity model. A natural extension is to consider a 
multisite impurity model ("cluster"), whose various self 
energies could be used to obtain a better representation of 
the lattice self energy. Several proposals have been made 
including a self-consistent embedding of a physical clus- 
ter ( "CDMFT" }i£ and a momentum space approxima- 
tion ("DCA")fii Recently a unifying "fictive impurity" 
( "FI" ) picture was presented , 12 ' 13 in which the different 
approaches were seen to correspond to different choices 
of basis in the same general expansion for the self energy. 

The relative merits of the different approaches have 
been debated?^ but, there have been relatively few com- 
parisions of the different methods in relevant physical 
limits. In this paper we take a step towards remedy- 
ing this deficiency by presenting, for the two dimensional 
half- filled Hubbard model in the strong correlation limit, 
a numerical study of real-space and momentum-space 
cluster dynamical mean field algorithms along with a 
comparison to analytics. A new feature of our analysis is 
that we are able to identify the contributions which arise 
from true intersite correlations (i.e. those not occurring 
in the single-site approximation) and compare them to 
exact (high temperature series) results, thereby quanti- 
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fying the degree to which the different methods capture 
the intersite correlations. 

Our results reveal that none of the methods give a par- 
ticularly good treatment of the intersite correlations The 
real space method discussed in Ref. [l2| has severe inad- 
equacies, which arise mathematically from an incorrect 
treatment of the pole structure of the self energy. The 
importance of respecting the pole structure of the self en- 
ergy was recently stressed by Stanescu and Kotliar. 14 Our 
results also point to a fundamental deficiency of the "Ac- 
tive impurity model" approach (in any of its implementa- 
tions): while general argument a 12 ' 13 guarantee that some 
cluster model exists which reproduces any given approxi- 
mation to a lattice model, the construction of the cluster 
model (in particular the choice of interaction terms) is 
not a trivial issue. While the DCA approximation pro- 
vides a better approximation to the intersite contribu- 
tions than do the other methods, none of the approaches 
are particularly accurate. We suggest that an impurity 
model with additional interaction terms would likely bo 
superior. 

The rest of this paper is organized as follows: section 
II defines the formalism we use and presents a few re- 
marks on issues related to numerical implementations. 
Section III presents our numerical results. Section IV 
gives analytical arguments which shed light on some of 
the findings and section V is a conclusion. 



which is the sum of a system-specific part and a univer- 
sal part, is minimized at the physical density and from 
which the ground state energy can be calculated. Den- 
sity functional theory became a useful tool following the 
demonstration of Kohn and Sharn^ that uncontrolled 
but reasonably accurate approximations to the univer- 
sal function could be constructed, and that a relatively 
convenient procedure for performing the minimization 
could be found. Similarly new progress in many-body 
physics has become possible following the formulation of 
an uncontrolled but reasonably accurate approximation 
to flskei along with a procedure for performing the min- 
imization. The approximation uj) — > S(w) (analo- 
gous to the local density approximation) was shown^ to 
permit the calculation of f2 s fcei in terms of the solution of 
a "quantum impurity model" with parameters fixed by 
the stationarity condition, Eq. |2]). 

The possibility of extending the approach to capture 
some part of the momentum dependence was alluded 
to in early work. 4 A discussion was given in previous 
work by some of usi 2 - (see also closely related work of 
Potthofft 3 ). In this paper we p resent detailed studies 
using the formulation of Ref. 12. To establish the nota- 
tion and define clearly the assumptions made, we outline 
the results of Ref. Qjl here. First, one approximates the 
momentum-dependence of the self energy in terms of a 
finite number, N, of basis functions 4>j(p) 



II. FORMALISM 



A. General Aspects 



A general result of many-body theory is that all elec- 
tronic physics of a given system can be obtained from the 
"Luttinger-Ward functional"— $ of the electronic self en- 
ergy E(p,w): 



$ = f7 sfcei ( ( T)-Trln[g - 1 -E] 



(1) 



Here Go — {dt — -Ho) -1 is the Green function of the as- 
sociated noninteracting model and the Luttinger-Ward 
functional Q s kei is defined as the sum of all vacuum to 
vacuum skeleton diagrams (with appropriate symmetry 
factors) and is here viewed as a functional of the elec- 
tronic self energy. The physical self energy correspond- 
ing to a given noninteracting problem (specified by Go) 
is determined from the stationarity condition 
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which follows because il s kei has the property that 



(2) 
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The situation is closely analogous to that obtaining in 
density functional theory, where general theorems 1 ^ guar- 
antee the existence of a functional of the electron density, 



E(p,cj) 
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such that as N — > oo T, approx — > E(p, uj). If one substi- 
tutes Eq. (0| into Eq. |T]) one obtains a functional ^approx 
of N self energy functionals Sj. The stationarity condi- 
tion Eq. ([2]) becomes the dynamical mean field self con- 
sistency condition 
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The most general such functional $ 



approx 



is an A-site 

quantum impurity model, which should be regarded sim- 
ply as a machine for computing the A functions £j(w) 
needed to generate the approximation for S(p, u>). The 
impurity model need not be a physical subcluster of the 
original lattice and is therefore referred to as "fictive" . 

Specifying the impurity model is not a trivial issue. 
The usual procedure is to assume that it is given by the 
action 



Simp = J dTdr'aij(T - T / )V>j(7" / )V ; i( r ) + 



Hi, 



(6) 



where Hi nt is exactly the interaction terms of the original 
lattice and the ay are mean field functions to be deter- 
mined from the self consistency equation. The impurity 
model is then some sort of self-consistently embedded 
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sub-cluster of the lattice model. Interactions extending 
outside the cluster are neglected. 

Reference [TH showed that the different cluster dynam- 
ical mean field schemes proposed in the literature are all 
variants of this general scheme, with the differences aris- 
ing from different choices of basis function <pj (p) . How- 
ever, while it is clear that as N — * oo the procedure 
converges to the full solution of the lattice problem, it 
is not clear that at any finite N the impurity model 
ansatz Eq. ([6]) generates the functional <E> which would 
be obtained by replacing E(p, u>) by T, approx {p, w) in fl ske i 
above. As we discuss in more detail in the conclusions, 
one interpretation of the results we present is precisely 
that the ansatz Eq. (jHJ) is not adequate. 



B. Models and Approximations 
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FIG. 1: Upper panel: Brillouin zone partition correspond- 
ing to 4-site dynamical cluster approximation. The zone is 
partitioned into four tiles, tile centered at momentum (0, 0) 
(unshaded), one centered at momentum (dark shading) 

and two at momenta (0,7r) and (w, 0). Use has been made 
of invariance under translations by integer multiples of 2n. 
Lower Panel: real space structure of corresponding cluster 
model 

For our specific computations we study the Hubbard 
model with nearest neighbor hopping on a two dimen- 
sional square lattice, make two and four site approxima- 
tions and consider three choices of basis function <b„. The 



first is the Dynamical Cluster Approximation (DC A), in- 
troduced by Jarrel and co-workers.— In present language 
the DCA corresponds to partitioning the Brillouin zone 
into a finite number of regular tiles (square, for the two 
dimensional square lattice we consider here) labelled by 
their central momentum Pj and choosing the basis func- 
tions (j)j(p) to be equal to unity of p is within the tile 
centered on Pj and to be zero otherwise. The partition- 
ing for the 4-site approximation is shown as the upper 
panel in Fig. [1] The result is a piecewise constant lat- 
tice self energy specified by the functions Es (w) giving 
the value of the self-energy in each Brillouin zone region, 
thus: 



(7) 



The corresponding impurity model is the four site cluster 
shown in the lower panel of Fig. [TJ This cluster has four 
self energies, corresponding to the on-site, first neighbor, 
and second-neighbor separations; these are related to the 



Ss via 

j 

E((0,0),o;) = So(a;)+2Si(w) + E 2 (w) J 
E((7r,0),w)=S((0,7r),w)=S (w)-E3H, (8) 
E((7r,7r),w) = E (w)-2E 1 (w) + 5i(w). 

The second choice of basis function is the 'CDMFT' 
approach introduced by Kotliar and co-workers In this 
approach one partitions the real space lattice into a pe- 
riod array of regular placquettes (supercells) as shown 
in cf. Fig. [21 so that the Hamiltonian becomes H = 
H p iac + T with H p i ac = + H int an impurity model 
defined by the hoppings and interactions on the plac- 
quette and T the interplacquette hopping. The cluster 
is treated as an impurity and is solved, leading to a self 
energy E which is a matrix in the space defined by the 
cluster. The lattice Green function is G _1 = cj — E(p) — E 
with E(p) = #o +T(p). 

The CDMFT approximation necessarily breaks some 
of the lattice symmetries. In the 2-site approximation 
both point group and translational symmetries are bro- 
ken. Various choices are possible. For the choice dis- 
played in Fig.[2]in which the unit cell is chosen so that the 
primitive translation vectors are u = x+y and v = x — y. 
Indexing u by i and v by j we have 
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(9) 
(10) 



while the interplacquette hopping connects, say, site 2 on 
placquette (i, j) to sites 1 on placquettes {i, j+l) 

and (i + 1, j + I), so that after Fourier transformation 
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0+1, j) 



would lead to a 2d + 1 site impurity model which could 
be solved to specify the quantities So, Ej, . . . However, if 
the point group symmetry is unbroken then many of the 
self energies are equal and one should be able to obtain 
the self energies from a smaller impurity model (two site, 
if only nearest neighbor terms are retained; four site if 
first and second neighor terms are retained). Defining, 
for the two dimensional case, 



7 P = ^ {cos(p x ) + cos (p y )}, 

,( 2 ) - 



jf J = cos(p x )cos(p y ). 
We write the self-energy for the 2-site cluster as 

E(p,w)=E +4o i> Ei(w), 
and for the 4-site cluster as 

E(p, u) = E + 4 7 «E!H + 4 7 ( 2 >E 2 H. 



(14) 
(15) 

(16) 

(17) 



In the four site case the cluster model to be solved again 
has the topology shown in the lower panel of Fig. [TJ 



C. Numerical techniques 



FIG. 2: Panel a): Possible partitioning of 2d square lattice ap- 
propriate for 2-site CDMFT. Panel b): Possible partitioning 
of 2d square lattice appropriate for 4-site CDMFF 



with V(p) = 1 + e iV2p " + + e <v^(p«+P«) anc j Puj Pv 

vectors perpendicular to v and u respectively. In the four 
site CDMFT method, 
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(12) 



while the interplacquette hopping connects, say, site 1 
on placquette to site 2 on placquette (i,j — 1) and 
to site 4 on placquette (i — l,j), so that after Fourier 
transformation 



E(P) 



1 





e -2ip y 





1 



„2ip y 



?-2ip* 







1 






. e 2i Py 



1 





e ~2ip y 





(13)" 

A third choice of basis function arises from a more 
straightforward real space expansion: E(R, ui) — > E(R = 
0, oj) + J2 a E(a, w) + . . . with a the set of vectors connect- 
ing a site to its neighbors. For historical reasons we refer 
to this as the (real space) FI method. Retaining only a 
few terms in this sum leads to a momentum space self 
energy expanded in the standard orthogonal harmonic 
functions. For example, if only on-site and nearest neigh- 
bor terms are retained, then a d-dimensional cubic lattice 



We used two numerical techniques to solve the TV- 
quantum impurity problem: quantum Monte-Carlo 
(QMC) using the Hirsch-Fye algorith m 4 ' 20 ! 21 and a re- 
cently formulated^ 2 - semiclassical approximation. 

The QMC technique is standard, but one technical is- 
sue requires comment. This method is formulated in 
imaginary time, and involves discretization so that the 
imaginary time integrations in Eq. ^ are approximated 
as the sum over the L 'time-slices' r n — n(3/L. The 
computation time scales as L 3 so the number of time 
slices which can be taken is limited and at lower tem- 
peratures (larger (3) the time step Ar = j3/L becomes 
uncomfortably large. The self-consistency step requires 
frequency-space information, and hence a Fourier trans- 
form which becomes inaccurate above the 'Nyquist fre- 
quency' ljm = 7rL//3. An additional difficulty is that the 
Green function has magnitude and derivative discontinu- 
ities across r = (corresponding to power-law decay at 
high frequencies).; these must be represented accurately 
to obtain the high frequency behavior correctly. Doing 
so is difficult because for in the strong interaction limit 
G varies rapidly near Ar = 0. Thus, the errors at fre- 
quencies of the order of the Nyquist frequency are large 
and for the range of L accessible to us the resulting errors 
are too large to yield reasonable estimates of the Green 
function. 

To mitigate the problems one must incorporate a priori 
information about the short time behavior of the Green 
function into the analysis, by using the short-time ex- 
pansion of the equation of motion for the lattice Green 
function to fix the size of the magnitude and derivative 
discontinuities across r = 0. This is typically done via 
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the following trick 2 ^: one introduces a "model function" 
G m odei{ T — T ') which has the correct high frequency be- 
havior up to some order uj~ m and considers the differ- 
ence 5G(t) between the model function and the QMC 
data. The low frequency behavior of the model function 
is not important; we took the appropriate momentum in- 
tegrals of the lattice Green's function with the self-energy 
E CT (w) = U[n- a - 0.5) + U 2 n_ a {l - n- a )/oj. The dif- 
ference function 5G is by assumption smooth near t = 
and in particular the first m — 1 derivatives are contin- 
uous. In practice a reasonable choice of model function 
leads to a SG which varies much less rapidly near r = 
than the original data or the actual Green function. 

By taking the difference between the QMC data and 
the model function, one obtains an approximation to 5G 
at the discrete points r n = n/3/L. One includes the a- 
priori information concerning the high-frequency behav- 
ior by performing an order m spline fit assuming that 
across r = the first m — 1 derivatives are continuous. 
In the single-site DMFT a cubic spline was found to be 
sufficient^ but in our investigations of multisite models 
it was found necessary to fix the a->~ 4 behavior of the 
Green function in order to control the high frequency be- 
havior of the first neighbor self energy. This necessitated 
the use of a fourth order spline fit to the QMC data. 
Fig- El demonstrates this effect, comparing the results of 
three different computations of the on-site self energy us- 
ing a two-site cluster (in the real space formulation) to 
the known high frequency behavior. 
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FIG. 3: Two-site cluster approximation to frequency depen- 
dent on-site self energy for paramagnetic phase of two di- 
mensional Hubbard model with U/t = 16 at temperature 
T/t = 0.2 calculated by standard procedure (third order 
spline fit; dash-dot line); four order spline fit directly to QMC 
data (dashed line); fourth order spline fit plus model func- 
tion subtraction (solid line). (The model function was ob- 
tained by appropriate integral of lattice Green function with 



T, a (uj) = U{n- a - 0.5) + U 2 n- a {l 



r)/w.) Results are 



compared to exact leading analytical high frequency result 
(dotted line). 

The QMC method remains very computationally ex- 



pensive; one requires a time slice short enough that 
UAt < 1 and very good statistical accuracy in the com- 
puted G's. To access a wider range of parameters we 
also used a semiclassical approximation we have recently 
developed which is much less computationally expensive. 
The SCA method is described in detail elsewhere,— so 
we mention here only a few aspects relevant to its imple- 
mentation in the present case. 

For an iV-site impurity model the partition function is 
defined as a functional integral over the 27V-component 
spin and site-dependent spinor fields and c as 



where 
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(19) 



with a the 2N x 2N matrix mean field function. To 
derive the semiclassical approximation we rewrite the in- 
teraction term as 



Vn jA (T)n u (T) = j[N*(T) 



Mf(r)] , (20) 



with n m = i ((n T + n x ) 2 - (n T - n x ) 2 ) = \{N 2 ~M 2 ). 
N is the number of particles and M is the magnetisation 
on the site. We then make the usual continuous Hubbard- 
Stratonovich transformation to decouple the M terms via 
a site-dependent auxiliary field <pj (r) which we assemble 

into an iV-component vector 4>. The semiclassical approx- 
imation is to retain only the zero-matsubara frequency 
term in the functional integral over </>. To this level of 
approximation the N term may be ignored because we 
work at half filling in a particle-hole symmetric model. 
We may then integrate out the electrons and obtain 



Z = J d$e s °" [a ' 4 



(21) 



where the effective action S e f f — (3V is defined by 

N , 



V(ct>) = ^\<j>\ 2 



T 2J Trln -a^{uj n ) - !<(> ■ a 



, (22) 



with 1 the 2N x 2N unit matrix. 

The integral over is a simple classical integral which 
may be done without too much difficulty. However, at 
strong coupling and low temperatures V is characterized 
by several very deep minima with high barriers between 
them and it convenient to make a further simplification 
and approximate the integration over <fi by the sum over 
the minima: 



N min 
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3=1 



-0V(4>j) 



(23) 
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where N m i n is the number of minima in potential V(<j)). 
This approximation corresponds to approximating the 
spins as Ising variables. 
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FIG. 4: Single-impurity DMFT results for the spectral func- 
tion Ao = — — ImGo computed by QMC and SCA. Solid line is 
QMC and broken line is SCA result. U/t = 20 and T/t = 0.5. 



The semiclassical approximation is reasonably good 
in the strong coupling regime. It reproduces all of the 
qualitative features found in the QMC calculations, and 
is reasonably quantitatively accurate. As an example, 
Fig. 2] shows the density of states calculated by analyti- 
cal continuation of QMC and semiclassical data for the 
single-inpurity Hubbard model. One sees that the semi- 
classical method places the Hubbard bands very close to 
the correct positions. Similarly, Fig. [5] shows the on-site 
and first neighbor spectral functions computed using the 
real space (upper panel) and DCA (lower panel) two-site 
approximation to the square lattice Hubbard model for 
the same parameters. Note that the unphysical feature 
in the density of states near u> = (to be discussed in 
more detail below) evident in the real space calculations 
but not in DCA is reproduced (or not reproduced) by the 
semiclassical approximation as appropriate, although the 
magnitude is not accurately determined 

To summarize, the semiclassical and QMC methods 
yield very similar results for the parameters relevant to 
this study. The semiclassical method is orders of magni- 
tude less computationally expensive. For example, per- 
forming one two-site cluster calculation at U/t = 20 and 
T/t = 0.5 required about 24 hours on a 2.4 GHz Pentium 
computer, essentially because the partitioned phase space 
means that up to 10 7 configurations must be generated 
to sample the entire phase space adequately. By contrast 
the semiclassical calculation requires about 5 minutes on 
the same computer. Therefore most of the results pre- 
sented below are obtained from SCA calculations. 




FIG. 5: 2-site Fictive impurity(upper panel) and 2-site 
DCA (lower panel) results for the spectral functions Ao = 
-ilmGo and Ai = -^ImGi computed by QMC and SCA. 
U/t = 20 and T/t = 0.5, paramagnetic order. 



III. NUMERICAL RESULTS 

A. Overview 

In this section we present numerical results obtained 
by the methods described in the previous sections, and 
we compare these to high temperature series results. We 
study four quantities: the local density of states 



N(w) 



1 f d 2 p 



7T J (27r) 2 



ImG(p, w), 



(24) 



the internal energy, given in the paramagnetic state by 
(the 2 is for the spin sum) 



E = 2 / — / ^4W( w )Im 



7T J (2tt) 5 



(25) 

the impurity model nearest neighbor spin-spin correlation 
function (a±a2) and the phase diagram. 
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B. Density of States 



C. Internal Energy 




FIG. 6: Spectral functions, obtained by FI method, DCA and 
CDMFT at U/t = 16 and T/t = 0.3, for l-(upper panel), 2- 
(middle panel) and 4-site(lower panel) cluster, with antiferro- 
magnetism suppressed so model is in the paramgnetic phase. 



Figure [6] shows the single particle density of states, cal- 
culated by maximum entropy analytical continuation of 
our numerical solution of the dynamical mean field equa- 
tions, for the paramagnetic phase of the square lattice 
Hubbard model with U/t = 16 and T/t = 0.3. (For 
the real-space approximation scheme this temperature is 
below the actual Neel temperature. In the data shown 
in Fig. [5] the magnetism has been suppressed to present 
results in the paramagnetic phase for all cases consid- 
ered) . The upper panel shows the spectral function com- 
puted from the single site DMFT; the model is obviously 
in the Mott insulating phase, with well separated upper 
and lower Hubbard bands. The middle panel shows the 
real space, DCA and CDMFT results for the density of 
states obtained from a two-site cluster. As in the single- 
impurity model, one observes the two Hubbard bands, 
the narrowing of the bands relative to the single- impurity 
case is a consequence of intersite magnetic correlations; 
indeed even in the one-site model, in the fully ordered an- 
tiferromagnetic case the bands are substantially narrower 
than in the paramagnetic phase. One also sees that in 
the real space (FI) method a small band of mid gap states 
exists. The lowest panel shows results obtained from four 
site clusters. One sees clearly in the FI method that the 
area of the mid-gap states decreases as the cluster size 
increases, and the frequency dependence changes. In a 
Mott insulator, the on-site self energy diverges as ui — » 0. 
The mid-gap states imply that in the FI cluster approx- 
imation S becomes small for some ui ~ 0. These results 
suggest that convergence to the infinite cluster size limit 
is not uniform in frequency. 



In this subsection we present results for the internal 
energy E = (H) computed from Eq. (|2"5")) . We remove 
the Hartree shift U/4 and the chemical potential. We 
compare the calculated results to analytical large U re- 
sults, which have been obtained up to 0(t 4 /U 3 )M^. To 
order t 2 /U one has 



t z 



£(2) =-2-tanh — 
U \ 4T 



U 



* 2 {^tanh(^ 



4T/ 



3} 



2Tcosh 2 



(26) 



E^> is shown as the light dotted line in Figure[7J It in- 
cludes terms from virtual excursions of an electron from 
one site to neighboring sites, but these average incoher- 
ently over the different relative spin orientations, so do 
not involve intersite correlations. 

In this model at half filling, the nontrivial inter- 
site physics is spin correlations and appears first at 
0(t 4 /U 2 T) ~ J 2 /T. To obtain results to this order we 
computed E = fl — TdVl/dT numerically from the ex- 
pressions for the thermodynamic potential Q presented 
by Kubo^ The result is plotted as a heavy dashed line 
in each panel of Figure [JJ 

Internal energy results as a function of temperature 
at U/t — 16 are shown in Figure [JJ for the real space 
(upper panel), DCA (middle panel) and CDMFT (lower 
panel) schemes, along with analytical results. For the dy- 
namical mean field methods, we show results both in the 
paramagnetic state and the antiferromagnetic state. The 
calculated Neel temperature is visible as the point of dis- 
continuity in the E(T) curves; for T < T/v we show both 
the antiferromagnetic state energy (lower curve) and the 
energy of the paramagnetic state (obtained by artificially 
suppressing the Neel state). We note that in order to ob- 
tain accurate energies the high frequency behavior of the 
Green functions must be carefully controlled. 

All of the curves display three temperature regimes: 
a very high-T regime (for the parameters considered 
here, beginning at T/t > 0.75 ) where the energy in- 
creases with increasing T, an intermediate T regime (here 
~ 0.5 < T/t < 0.75) where the energy is approximately 
T- independent, and a low-T regime in which the energy 
exhibits a strong T dependence. The increase of E with 
T in the high-T regime arises from real thermal excita- 
tions over the Mott-Hubbard gap [cf. the second term in 
Eq. (f2"fJ]) ]. The more rapid upturn of the DMFT results 
relative to the series expansion is an artifact of the SCA, 
which overestimates the effect of thermal flucuations on 
the gap. This regime will not be discussed further here. 

In the intermediate T regime, the excitations into the 
upper Hubbard band are quenched, and intersite spin 
correlations are slowly developing. The single-site DMFT 
neglects intersite correlations entirely in the paramag- 
netic phase; thus in this regime the single-site DMFT 
result is essentially independent of temperature, and is 
seen to be very close to the second order series result 
t/U = 0.0625; we would expect corrections to be of rela- 
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FIG. 7: Internal energy E/t as a function of temperature ob- 
tained by FI method (upper panel) , DCA (middle panel) and 
CDMFT (lower panel) at U/t — 16 and compared to analyti- 
cal results. Fat solid curves are for single impurity, solid(solid 
with squares) curves for 2-site(4-site) cluster in antiferromag- 
netic state, dashed(dashed with squares) curves for 2-site(4- 
site) cluster in paramagnetic state, stars and crosses are ob- 
tained as described in the text from the large U expansion, 
to the order in U indicated. The rapid rise with temperature 
of the DMFT results for temperatures T > t is an artifact of 
the implementation of the semiclassical approximation based 
on Eq. (|23p used here, which overestimates the contribution 
of excitations across the upper Hubbard band. 



tive order 1/U 2 ~ 10 -2 , essentially invisible. The effect 
of intersite correlations is visible below the Neel temper- 
ature. 

Both 2 and 4-FI methods produce energies which lie 
above the single-site DMFT curves and which increase 
at low T. We conclude that in these methods the in- 
tersite spin correlations are wrongly treated, leading to 
an 0(J 2 /T) contribution to the energy with the wrong 
sign. The physical origin of the error is the mid-gap 
states which shift weight of order £ 4 /<7 4 in ImG from the 
vicinity of the lower Hubbard band up to the chemical 
potential, thereby raising the energy. 

The CDMFT and DCA methods produce energies 
which lie below the single-site curve, indicating that they 
provide a qualitatively correct treatment of the inter- 
site spin correlations. The quantitative accuracy may 
be judged from the separation between the DMFT cal- 
culations and the series results. The agreement is not im- 
pressive. The CDMFT intersite energy is far too small, 
while in the DCA method the 4-site cluster produces an 
energy in worse agreement with the correct answer than 
does the 2 site cluster. 

Finally, we note that in the four-site methods, if the 
Neel ordering is suppressed an apparent first order tran- 
sition (most probably to a dimerized spin state) occurs. 



D. Neel temperature 

The Neel transition temperature Tn was identified 
with the temperature corresponding to the kink in the 
antiferromagnetic E(T) curve. We note that our meth- 
ods are mean field methods. We have verified the values 
by writing an independent code to obtain the tempera- 
ture dependence of the staggered magnetization. In the 
two dimensional models studied here spatial fluctuations 
drive Tn logarithmically to zero; for the small clusters 
studied here our computed Neel temperature is therefore 
best interpreted as a scale below which the spin-spin cor- 
relations become appreciable. The computed mean-field 
phase diagram is shown in Fig. [5] In the small U limit 
all the methods agree reasonably well with each other 
and with the simple analytical results. This finding is 
in agreement with a detailed study of the size depen- 
dence of the Neel temperature at small J7i2& However, 
at large U substantial variation exists. We observe that 
the single-impurity calculation produces results in much 
better agreement with the large-U limit than the cluster 
methods. The FI method grossly overestimates Tjv. We 
believe that the overestimate occurs because the order- 
ing eliminates the mid-gap states, thereby substantially 
lowering the energy, (cf. Fig. [7]). The unphysical nature 
of the FI results means that computations of the four- 
site FI method are not worth performing and are not 
presented here. 
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FIG. 8: Neel temperature Tjv/t vs. on-site interaction U/t. 
Solid curves are fits of the data. 



E. Impurity model spin correlations 

We finally consider the spin correlations in the impu- 
rity model. (Note that the "Active" nature of the im- 
purity model means that it is not to be thought of as a 
physical subcluster of the lattice, so the relation of the 
impurity model spin correlations to the actual spin cor- 
relations in the lattice is not entirely straightforward.). 
In Fig. [5] we show the comparison of the NN spin-spin 
correlation to the 2- and 4-site CDMFT (lower panel), 
DCA (middle panel) and FI method (upper panel) re- 
sults as a function of temperature. Also included is the 
leading term (oxCTa) = —t 2 /(TU) in the appropriate high- 
temperature-series expansion. We see that the various 
methods obtain results which have the correct tempera- 
ture dependence, but with magnitudes somewhat at vari- 
ance with the exact results. We observe that the un- 
derestimate of the intersite contribution to the energy 
is not reflected in an underestimate of the cluster spin- 
spin correlations, suggesting that the deficiencies of the 
methods have to do with interactions which extend out- 
side the cluster considered. We also note that for the 
sizes available to us, increasing cluster size does not lead 
to improved agreement. 



IV. APPROXIMATE ANALYTICAL 
TREATMENT 

A. General formulation 

In this section we present approximate analytical cal- 
culations which provide some insight into the numerical 
results. The calculations are based on an approximation 
to the semiclassical method of Ref. |22j This first sub- 
section gives some general considerations. The next sub- 
section presents the relevant aspects of the approximate 
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FIG. 9: Nearest neighbor spin-spin correlation as a function 
of temperature obtained by FI method (upper panel), DCA 
(middle panel) and CDMFT (lower panel) at U/t = 16 (see 
also section WV} . Also shown is the high temperature series 
result for the Ising approximation to the square lattice Heisen- 
berg model (<j\ z (T2z) = —t 2 /(TU), which is the result to which 
it is appropriate to compare the numerical calculations per- 
formed using Eg 1231 
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solution of the impurity model (which is the same for all 
methods). Subsequent sections combine these formulae 
with appropriate self-consistency conditions to obtain re- 
sults for the single-site model and the 2 and 4-site DCA, 
CDMFT and FI approaches. 

In developing the analytical approximations it is useful 
to alternate between two basis choices for the impurity 
model. First, a real-space basis with on-site ao and inter- 
site dj^o mean field parameters. The key simplification 
of the large U half-filled limit is easily seen in this basis: 
the magnitude of the intersite terms a^o is much less 
than | dp — (f> 2 1 . Assuming no spatial symmetry breaking 
(so 1 0| is the same on all sites) and expanding to second 
order in the parameter af^ /(a,Q — 4> 2 ) leads to an expres- 
sion involving the mean field parameters and the intersite 
spin correlations, which may be treated analytically. 

For solving the self consistency condition it is more 
convenient to consider the impurity model eigenbasis. An 
AT-site impurity model involves an N x N matrix mean 
field function a, Green function G and self energy X re- 
lated by £ = a G _1 . We restrict attention to the 
paramagnetic phase so a, G, S are proportional to the 
unit matrix in spin space. For the models of interest the 
orbital-space matrices may be simultaneously diagonal- 
ized, so that for the N eigenmodes A we have 



G\ = (a\ - Da) 1 



(27) 



B. Impurity Model and Self-consistency condition: 
large U limit 

An N site impurity model is specified by a set of P + 1 
N x N matrices Mj. The impurity model action Si rnp is 

p 



Simp — ^ ^ 



(29) 



For all models, M is the N x N unit matrix 1. For 
the two-site model, Mi = t x and the eigenvectors are 
correspondingly the even and odd combinations 



1 / x 
a e = -j={a + ai), 

a Q = —7= (ao - ax) . 

For the four-site model, P = 2 with 



Mi = 




M 2 



The eigenvectors are 



10 
1 
10 
10 



(30) 
(31) 



(32) 



(33) 



The DMFT self consistency equation are obtained by 
relating the lattice and impurity model Green functions. 
Different schemes involve different methods for relating 
the impurity model Green function and self energy to the 
lattice green function and self energy. In the impurity 
model eigenbasis we have 



\S) = 



\Y) 




\X) = 



\D) 




ax 



Sa = 



(dk) 



1 



iuj — Ek — to) 



(28) 



Here the notation [/ (dk)\ . denotes the details required 
for the particular DMFT scheme. We may then ex- 
pand the right hand side, noting that in a Mott insulator 



\iu) - S| > \e k \ and that E(fc,w) = S A - 

(1.2) 



„(1) _ U 2 , „(2) 



with s y ^'^' small. This formulation enables one to solve 
for the mean field parameters without explicitly comput- 
ing the Green functions or the sub-leading contributions 
to the self energy. 

In the rest of this section we present the details of the 
large-U analysis. We first give the analytical solution of 
the general impurity model, then present the connection 
to the lattice, and finally give results for physical quan- 
tities. 



so 



as 
ax 

fly 

an 



a-i, 



a + \pia\ 
a - a 2 , 
a - a 2 , 
a — V~2ai + a 2 



(34) 
(35) 
(36) 
(37) 



To solve the impurity model we proceed from Eq. (|23j) 
in the large U limit. We treat the integral over the mag- 
nitude of the auxiliary field by the steepest descent ap- 
proximation, so that in the large U limit \(f>\ w U/2 is the 
same on each site but the direction flj may vary. The 
partition function is then 



(dQj)e' 



-/3Ve//({%}) 



(38) 



with the \4>\ fixed by 



ov 
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We consider temperatures low enough that thermal 
excitation into the upper Hubbard band may be ne- 
glected (mathematically this means we replace Tj2^ n 
by J duj / '(2tt) in all expressions, with one exception dis- 
cussed below). 

Hubbard-Stratonovich transformation followed by in- 
tegration over the fermion fields leads, in the semiclassi- 
cal approximation, to 



For comparison to the numerics we note that in the 
Ising approximation used in the numerical calculations 
the factor 3 in the denominator of the right hand side of 
EqgHis absent. 

The impurity-model Green functions are Qj"- p = 
Qj'M.j with gj = S In Zj (2NSa,j). In both two and four 
site models we find (assuming a 2 = 0) 



S;. 



Tr In 



E a i M i 

j=0,P 



N<t> 2 

IT 



(39) 



with D a diagonal matrix with entries (f>fli -a. At large U 
and half filling we have | a§ — 2 1 ^ |ojyo| 2 - Expanding 
to second order gives 



Si; 



N<t> 2 



iVTrln [al - <j> 2 } 



- E ; T^^2 Tr t M * to 1 + D ) M ^ M + D )l -( 4 °) 

Taking the trace explicitly yields, for two and four-site 
models, 



S 



„„„ -^- + 2Trln[a 2 -^ 2 ] 
#2 1 4Trln[a 2 -<£ 2 ] 



2a 2 [a 2 



(a 2 - <f) 



2 fii -tl 2 ' 

—s—&) 



imp 

a 



U 



2 {4a 2 + 2 ( 



f2l • fl 2 + ^2 • ^3 + ^3 ' ^4 + ^4 ' ^1 



)} 



|2a 2 + 2<j> 2 {ti x -n 3 + n 2 - f2 4 ) I 



We shall see that for the Hubbard model with nearest- 
neighbor hopping, a 2 = to the order to which we work. 
In this case, for both two and four-site models, the mean 
field equation fixing <j) is 



— -^E 



1 a 2 {2a 2 + S ' 



<t> 2 )) 



a 2 -^ 



with S the nearest neighbor spin correlation given for 
N = 2, 4 by 



S=(Sli-fk 



--E 



n{*l-4> 2 Y 



(44) 



where the second approximate equality comes from ex- 
panding Z to leading order in a 2 /(ag — 4> 2 ) and applies for 
T sufficiently greater than J = t 2 /U. Note that Eq|44lis 
written for Heisenberg spins; the semiclassical numerical 
method used here amounts to an Ising approximation in 
which the prefactor becomes A/N. 



,9o 



fJi 



ao L 4 («§ + 2 (l + 2^)) ' 
«o - <$> 2 \ {a 2 - <b 2 f 

-ai (ag + fS) 

K-0 2 ) 2 



(45) 



(46) 



General expressions for the self energy are cumber- 
some. By combining go and g± into the appropriate 
impurity-model eigencombinations we find that 



(47) 



In the low frequency limit \u>\ -C cf> we have, for the two 
site model 



a Q + axS' 

j3 



S = 



ao — aiS" 
while for the four-site model 

S s = 



ao + a\S 



(42) 



a 

E - 02 
- — , 

a 

V D = 



ao — a\S 



(48) 
(49) 

(50) 
(51) 
(52) 
(53) 



For the impurity models in the insulating regime, we 
will find that a ~ u> while ai ~ t. Thus the low fre- 
quency behavior of the impurity model self energies is 
well approximated by the simple pole 



(54) 



In the single-site dynamical mean field approximation, 
Q\ = but in general Q\ is of order t with a prefac- 
tor which depends on the intersite spin correlations and 
becomes very small at T > t 2 /U = J. 

Differences between dynamical mean field schemes 
arise from different ways of combining the impurity 
model self energies into an approximation to the lattice 
self energy. In the DCA and CDMFT approaches, the 
impurity model self energy translates essentially directly 
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into a lattice self energy, so that the pole structure is pre- 
served. In the FI approach, the situation is different. For 
example, in the two-site model one has, at low frequency, 



Ef/M « (f> 2 



l + 2d lk , l-2d 7 * 



(55) 



with Q e Q Q . Eq. (|55[1 implies that at a general k the ap- 
proximate self energy has two poles with a zero-crossing 
between them. This incorrect analytical structure leads 
to the midgap states found numerically 



C. Single-site approximation 

In the single-site problem, the on-site terms are the 
only ones present so we set a\ = S = in the formulae of 
the previous section. The impurity model Green function 
is 

a 



G, 



so that 



E = — . 

The self consistency equation is (dk) = d d k/(2ii) 

-l 



ao - S = 



(dfc)- 



l 



(56) 



(57) 



(58) 



Now in a Mott insulator we expect \iu> — £| 3> |£fc| 
Thus we rewrite Eq. ([55]) as 



«o 



(icj — E) 



1 



(iw - E) { 1 + 



A',, 



(iw — E) z 



where 



Thus 



(dk)e 2 k = 2dt 2 . 



a 
E 



-A2 



K d 



1 



LJ 2 + (j) 2 

K d 



(59) 

(60) 

(61) 
(62) 



uj 2 + 4> 2 / 

while substitution into Eqs. (|43l25p . expansion and re- 
placement of the frequency sums by integrals gives 

, = U_K± 
^ 2 2U' 

U K d _ U dt 2 
' ~8 ~ 4[7 ~ ~ ~8 ~ 2JJ' 

We observe that to this order in the t/U expansion the 
single-site DMFT is in agreement with the exact result, 
Eq. 



E 



(63) 
(64) 



D. DCA 

In the DCA one covers the Brillouin zone with N tiles, 
A, which correspond to the eigenvectors of the impurity 
model. The self consistency equations are 



(dk) - 



-i -1 



(65) 



where J*, (dk) denotes an integral over tile A of the Bril- 
louin zone, normalized so J x (dk) = 1 An analysis identi- 
cal to that leading to Eq. ([55)1 gives, up to corrections of 
order t 3 /U 2 



a x =iw- I\- 



K X -I 2 



with 



(dk)s k , 



K x = J(dk)e 



(66) 

(67) 
(68) 



Note that in the limit spatial dimensionality d — > oo 
K ~ d whereas ^ A ^ 1 so that in this limit the model 
reduces to single-site dynamical mean field theory. 

In the two-site DCA the two eigenstates are even (e) 
and odd (6) and we find (in d = 2) 



I, 

K,. 



r(2) 



(2) 



implying 



Kn = K 



iuj { 1 + 



it 



-1.62i, 



K m _ (^y 



UJ 2 + (f> 2 



m 

aj = -J e = —5-. 



In the 4-site DCA we have 



Is 
Ix 

K s 
K x 
Let us define 

j(4) = 



Q4. 

-In = re -2Mt, 

7T 



K 



D 



0, 



At 2 



m 2 



= k y = u 2 - 



m 2 



7.24t 2 



0.76t 2 



(69) 
(70) 



(71) 
(72) 

(73) 
(74) 

(75) 
(76) 



' (Is ~ Id) = — » 1.80t, (77) 



2^/2 

= 1 y K\ = 4i 2 . 

A ^ — ' 



(78) 
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Then 



ao = = iu> I 1 
ai = 1^. 



W 2 + 2 



(79) 
(80) 



Thus for the N = 2, 4 site models the Ising version of 
Eq. (gH) implies 



I 2 

' N<t>T 



2/2 
NTU 1 



(81) 



with J given by either /( 2 ) or as appropriate. From 
Eq. B3]) we have 



17 



00 



- 2 (l + ^) +0 2 



7 2 {-2cO 2 + 5(-^ 2 +0 2 )} 



(cC 2 + 02)3 



r°° du 


i 






2 ) 



1 - 2- 



^ + 2 r 

7 2 0{-2^ 2 + 5(-^ 2 + 2 )}" 



(c 2 + 2 ) 3 



r°° duj 


r * { 




c 2 + 2 1 


( 


-a, 2 + 2 )" 


(uj 2 


+ 2 ) 3 



1 - 2- 



i^ 2 



(c 2 + 0T 



(82) 



2 ~ 2C/ 2 + 2f7 : 

with K given by X^ 2,4 ^ as appropriate. 

Finally, we consider the energy. Within DCA we have 

e v + i£ A (a; n ) 
yvi~ — 

n,X Jlx 



E DCA = 2T^2 I (dp) 



iuj n — S p — EaK) 



2T^ 



1. 



■1 + ( iuj n - -Ex ) G\(iu} n ) 



(83) 



We now rearrange Eq. (|83jl into a form more con- 
venient for the strong coupling expansion. We write 
G\ = (a\ — £\) 1 and by adding and subtracting ob- 
tain 



e =ty: 



(icu n G Q - 1) + — (iuj n - a x ) Gx(iv n ) 



(84) 

Finally, we note that because the same change of ba- 
sis diagonalizes G and a and G = J2a=o G«M„ and 
similarly for a with Tr [M;Mj] = NSij we have 



e =ty: 



N-l 



{a Q G - 1) + 2 (iuj n - o ) G - ^ a h G b 



6=1 



(85) 



Here each of the three terms is convergent at large u and 
the second and third are explicitly of order t 2 /U. 

We consider the three terms in turn. The first term is, 
explicitly 



E<» = T^(aoGo-l) 



a 2 -cp 2 



= T E 



a 2 -0 2 

U K -I 2 
7 + 2U 



i | a 2 { a 2 + 2 (l + 25)} ' 

(a 2 -0 2 ) 2 
| a 2 a 2 {a 2 + 2 (1 + 25)} 
(a 2 -0 2 ) 3 



(86) 



Similarly, use of Eq. (fTTj) gives 

£( 2 ) = T^2( lWn - ao )G 

n 

-2cc 2 (K - I 2 ) 



(a 2 -0 2 r 



. (87) 



while 



I 2 



Thus the total energy for the N = 2, 4 site DCA approx- 
imation is 



N 

DCA 



u 


K 


(jW) 2 5 


4~ ~ 


2U" 


2f7 


u 


K 


( 7 W) 4 


7 ~ 


2U ~ 


NU 2 T ' 



so that 



E 



2-DCA 



E 



4— DCA 



— 3.45 
— 2.62 



t 4 



U 2 T 7 

t 4 
'C/ 2 T' 



(89) 

(90) 
(91) 



We see that both two and four site DCA approximations 
lead to an expression for the energy which reduces to the 
single site expression if the spin correlation 5 = 0. The 
differences between the two and four site approximations 
have a small contribution from the difference in the fac- 
tors I but this is overcompensated by the factor of N 
in Eq. (|4"4"]). Numerically the coefficient of the 1/T term 
is seen to be larger for the two-site DCA than for the 
four-site DCA, so that (in agreement with the numerical 
results) the four site DCA is seen to have a sightly worse 
intersite energy than the two site DCA. 
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E. CDMFT 



FI Model 



The CDMFT approximation may be treated in a man- 
ner very similar to the DCA. The lattice Green function 
is a matrix in the space of the cluster states, so the self 
consistency equation is 



G imp = J {dk) [iu - £ - E(fc)]" 1 , 



(92) 



where the prime denotes an integral over the reduced 
Brillouin zone appropriate to the real-space tiling and 
the dispersion matrix E was given above in Eqs. (jlll 
[T3| . Expanding and noting that \iuul — S| ^> E(p) and 
that S is diagonal to leading order in t/U we find 

G imp = (tuj-Xy 1 (l + 1 (iw - £) _1 + K (iw - £)(9§) 

with 



K 



(dfc)E(fc), 



(dk)E(k) 2 



(94) 



(95) 



Thus, inverting once more and using again that X is ap- 
proximately diagonal we obtain 



a = itul — I — 



K-I 2 



In the two-site CDMFT we have 

1 



I = -t 



K 



-At 



1 



2/10 



(96) 



(97) 



1 



while in the four-site CMDFT we have 

I = -V2iMi « -1.4iMi, (99) 



K = -At 1 [ Mo + -M 2 



The solution of the self consistency equations and the 
analysis of the energy goes through as before; the only 
difference is in the values of the intersite parameters a\. 
We find 



j-l-CDMFT _ _. 
j^CDMFT = 

so the intersite term in the energy is 

i 4 



£j£j2— CDMFT 
^j^i— CDMFT 



-0.5 



t 4 



U 2 T' 



U 2 T 



(101) 
(102) 



(103) 
(104) 



Thus the CMDFT method underestimates the intersite 
correlations by a larger factor than the DCA but moves 
in the correct direction with cluster size. 



The analysis of the FI equations is not quite as 
straightforward as was the analysis of the DCA and 
CDMFT equations. We specialize at the outset to the 
two-site problem, which reveals the essential difficulties. 
In this case, the lattice self energy is (for the nearest- 
neighbor hopping model studied here) 



T,(k,u) 



1 + 2d lk 



1 - 2d lk 



(105) 



and the self-consistency equations are (for general d) 



G, 



1 

e.o 



(dk)- 



l± 7fc 



T,(k,u) - e k 



(106) 



Unlike the previously considered cases, the self en- 
ergy has a momentum dependence which interacts with 
the momentum dependence arising from the dispersion. 
Because S e and E D have poles at different energies [cf. 
Eqs. (|48l49p ], T,(k,oj) generically has two poles (with fc- 
dependent strengths and (except at special fc-points) a 
zero crossing between them, this structure is physically 
incorrect (the self energy should have only one pole at a 
given k) and the concomitant zero crossing produces the 
mid-gap states. 

To analyse the equation, say, for G e we write 



S(fc, uj) = H e (ij) + 



(1 - 2d 7 fc) , 



(107) 



assume the second term is small compared to the first and 
proceed as before. We obtain (in spatial dimensionality 
d) 



iw^f t 



Thus 



K d (i - f: 



<2l = t, 



1 

2d, 



(108) 



(109) 



and (again for the hypercubic lattice with nearest neigh- 
bor hopping, and keeping only terms up to order t 2 /U 2 ) 



2dt 



( 1 + f) ! 



1 2d 



(110) 



In the derivation of the single-site DMFT equations the 
d — * oo limit is taken with dt 2 held constant. In this 
limit, S ~ t 2 /(TU) ~ 1/d vanishes and the equations 
revert to the usual single-site DMFT form. 

Equation (|110|) is valid for T > ^fji, but the solution 
changes character for u> < \Jt 2 U /T. At high frequencies 
we may solve iteratively, obtaining 



ao^ioj <1 + 2dt 2 1 



1 

2d 



(111) 
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Thus if S is sufficiently small we find ao ~ lo and a\ ~ 
t. At lower frequencies, the structure of the equation 
becomes more complicated, because of the presence of 
Si ~ 1/ojq on the right hand side of the equation. This 
behavior arises because of the inappropriate combination 
of poles and produces the midgap states discussed above. 



V. SUMMARY AND DISCUSSION 

In this paper we have examined several multisite exten- 
sions of the dynamical mean field method in the strong 
coupling limit, which has not been the subject of previous 
systematic study. We have computed a variety of physi- 
cal quantities and compared these to available and newly 
comoputed analytical results. We were able to isolate the 
contributions which arise from nontrivial intersite (in this 
case, spin-spin) correlations. We found that an incor- 
rect treatment of these in a real space (FI) scheme pro- 
duces unphysical mid-gap states in the Mott insulating 
phase, and thus wrongly estimates the internal energy, 
Neel temperature, and spin correlations. The 'DCA' and 
CDMFT schemes did not lead to mid gap states, and 
produced results which are qualitatively correct. How- 
ever, substantial quantitative differences exist between 
the CDMFT/DCA results and the exact answers. 

^From a mathematical point of view the central diffi- 
culty with the FI approach is the pole structure of the 
self energy function. The importance of the pole struc- 
ture was stressed by Santescu and Kotliar— . In a Mott 
insulator the equation co — s p — S(p, u>) = has no so- 
lutions at low 00; the lack of solutions arises because the 
lattice self energy has the form given in Eq. ([53)1 : a simple 
low-frequency pole at each p. All of the DMFT schemes 
involve approximating the lattice self energy S(p, lo) by 
a combination of the N self energies T,\(u>) of an N site 
impurity model. Each of the impurity model self energies 
exhibits a pole at some low frequency The FI method 
combines the impurity model poles in such a way that at 
typical fc-values the lattice self energy contains N poles 
with zero crossings between them. This structure leads 
to mid-gap states. The DCA and CDMFT methods, on 
the other hand, translates the cluster self energy directly 
to the lattice, leading to a piecewise constant self energy 
with only one pole at each k, and therefore to no mid-gap 
states. 

An approximate analytical examination of the equa- 
tions in the strong coupling limit provides some addi- 
tional physical insight into the multisite DMFt method. 
At temperatures low enough that real excitations across 
the Mott-Hubbard gap may be neglected, the expansion 
may be thought of as sampling virtual excursions of an 
electron, which starts from one site, samples some num- 
ber of near neighbors, and returns to its starting point. 
The result depends on the intersite spin correlations. 
We found that all methods reproduce exactly the lead- 
ing 0(t 2 /U) result for the internal energy, but both the 
multisite methods provided incorrect and indeed in some 



cases unphysical estimates of the 0(t 4 / (TU) 2 ) terms. 

The correct value of the 0(t 2 /U) term has an in- 
teresting implication. An early examination of pos- 
sible multisite extensions of the DMFT method by 
Schiller and Ingersent^l has been interpreted as show- 
ing that straightforward cluster methods (such as the FI 
method) are fundamentally flawed, because they neces- 
sarily double-count processes involving the hopping of an 
electron from one site to another. Our finding, that all of 
the cluster methods agree with exact results at 0(t 2 /U) 
and that the disagreements arise from terms involving 
intersite spin correlations, calls this interpretation into 
question. It is obvious from our results that the various 
methods have various levels of flaws, but it appears that 
a fundamental overcounting is not among them. Instead, 
the errors arise from an incorrect treatment of the terms 
physically arising from intersite spin correlations. 

Additional insight into this question is provided by a 
strong coupling expansion performed for the Hubbard 
model by Pairault, Senechal and Tremblay^S, These au- 
thors presented results for the electron Green function up 
to third order in t. Because our quantity S ~ t 2 , their 
result for Go is equivalent to our result for this quantity 
with S = 0; at this order the cluster results agree with 
the one-impurity result. However, the results of Pairault 
et al imply that 

uj 2 4- ^M- 



The FI method obtains Eq. (fTT2"]) but with the coeffi- 
cient 3/4 replaced by 1/2 while the 2-site DCA method 
replaces the pref actor t by Ij, = 1.6t in d = 2 and the co- 
efficient 3/4 by Ij/12t 2 « 0.65. The differences between 
the DMFT and exact results arise from an inaccurate 
treatment of intersite correlations in the DMFT. 

Reference [H also showed that the strong-coupling ex- 
pansion for the Green function was not uniformly conver- 
gent, but in order to yield finite results had to be carried 
to an order which increased arbitrarily as the frequency 
was lowered. Our results show something similar: the 
'FI' method does not converge uniformly to the exact 
result as a function of cluster size and frequency or tem- 
perature. In the present case we traced the difficulty to 
mid-gap states induced by an incorrect approximation to 
the pole structure of the self energy. The other DMFT 
methods discussed here lead to self energies with the cor- 
rect pole structure, but to values for the intersite contri- 
butions to the energy which are in poor agreement with 
exact analytical results. The methods may be thought 
of as arising from resummations of particular classes of 
terms in the strong coupling expansion. The poor agree- 
ment with exact results suggests that the resummation 
is not precisely correct and indeed not necessarily partic- 
ularly accurate. 

We note that the weak point of the general arguments 
establishing the multisite DMFT approach is the choice 
of interaction terms in the impurity model. These are 
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always taken to be the same interactions as in the lat- 
tice model. We speculate that in order for the impurity 
model to represent the Luttinger-Ward functional with 
the truncated self energy, it may be necessary to incorpo- 
rate additional interaction terms, representing the effects 
of otherwise neglected intersite processes. In the model 
studied here the intersite processes have to do with spin 
correlations. The incorrect values of the intersite energy 
go along with more reasonable estimates of the cluster 



spin-spin correlations. This suggests that the difficulty 
with the energy relates to effective interactions which ex- 
tend beyond the cluster considered. . 
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